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We describe a first-principles statistical mechanics approach enabling us to simulate the steady- 
state situation of heterogeneous catalysis. In a first step density-functional theory together with 
transition-state theory is employed to obtain the energetics of all relevant elementary processes. 
Subsequently the statistical mechanics problem is solved by the kinetic Monte Carlo method, which 
fully accounts for the correlations, fluctuations, and spatial distributions of the chemicals at the 
surface of the catalyst under steady-state conditions. Applying this approach to the catalytic oxi- 
dation of CO at RuO2(110), we determine the surface atomic structure and composition in reactive 
environments ranging from ultra-high vacuum (UHV) to technologically relevant conditions, i.e. up 
to pressures of several atmospheres and elevated temperatures. We also compute the CO2 formation 
rates (turnover frequencies). The results are in quantitative agreement with all existing experimen- 
tal data. We find that the high catalytic activity of this system is intimately connected with a 
disordered, dynamic surface "phase" with significant compositional fluctuations. In this active state 
the catalytic function results from a self-regulating interplay of several elementary processes. 

PACS numbers: 82.65.+r, 68.43.Bc, 82.20.Uv, 68.47.Gh 



I. INTRODUCTION 

Rational design and advancement in materials science 
will ultimately rely on an atomic-scale understanding of 
the targeted functionality. In macroscopic systems of 
technological relevance, this functionality typically re- 
sults from the interplay of a large number of distinct 
atomic-scale processes. This makes quantitative calcu- 
lations challenging, since it requires not only to com- 
pute many individual elementary processes with high ac- 
curacy to reach the aspired predictive power, but also 
to appropriately account for the statistical mechanics of 
the interplay of all these processes. Over the last years, 
modern multi-scale methods that suitably combine first- 
principles electronic structure calculations with thermo- 
dynamics and statistical mechanics have advanced suffi- 
ciently to make this goal feasible. 1 

In heterogeneous catalysis the targeted functionality is 
the efficient conversion of a supply of chemicals A and B 
at the surface of a solid catalyst material into a product 
C that is then transported away. Since this corresponds 
to a thermodynamic open system, even steady-state con- 
ditions, where the production of C proceeds at a stable 
rate, are determined by kinetics. A quantitative compu- 
tation of the steady-state rate, measured typically as a 
turnover frequency (TOF, number of product molecules 
per unit area and second), requires therefore to explic- 
itly follow the time evolution of a large enough surface 
area, fully treating the interplay of all relevant underly- 
ing atomic-scale processes. For the example studied in 
this paper, the atomistic processes include the motion 
of the gas-phase molecules, dissociation, adsorption, sur- 
face diffusion, surface chemical reactions and desorption. 
Since many of these processes are activated and thus rare, 
macroscopic time scales need to be simulated, to arrive 



at meaningful conclusions concerning the effect of the 
statistical interplay. Thus, although individual processes 
occur on picosecond time scales, this may mean in prac- 
tice to simulate the system over time scales up to the 
order of seconds or longer. 

An additional level of complexity is encountered when 
aspiring a predictive, quantitative modeling that is based 
on understanding, and that is applicable to realistic en- 
vironmental situations of varying temperatures and pres- 
sures. This then not only excludes the use of empirical or 
fitted parameters, in which often several not further spec- 
ified processes are effectively "lumped together" , but also 
excludes the use of a mean-field description as e.g. mi- 
crokinetic rate equation approaches^. As we will further 
qualify below, parameters with clear microscopic mean- 
ing and a full consideration of the correlations, fluctu- 
ations and spatial distributions of the chemicals at the 
catalyst surface are crucial to properly understand and 
potentially rationally design macroscopic catalytic func- 
tionality. Each individual elementary step must thus at 
first be treated separately, and then combined with all 
the others within an appropriate framework. 

In a recent Letter we presented a first-principles statis- 
tical mechanics setup suitable to tackle this challenge 3 , 
and in this paper we give a detailed account of the em- 
ployed methodology. We use density-functional theory 
(DFT) together with transition-state theory (TST) to 
accurately obtain the energetics of all relevant processes, 
and subsequently solve the statistical mechanics problem 
by kinetic Monte Carlo (kMC) simulations. This two- 
step approach enables us to gain microscopic insight into 
the system, following its full dynamics from picoseconds 
up to seconds. In particular, the approach is employed 
here to provide an ab initio description of the compo- 
sition and structure of the catalyst surface in reactive 
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environments ranging from ultra-high vacuum (UHV) to 
technologically relevant conditions with pressures of the 
order of atmospheres and room temperature or higher. 
We also quantitatively compute the catalyst activity in 
terms of TOFs under all of these conditions. 

In the earlier Letter we used the CO oxidation at 
RuO2(110) to illustrate the new quality and novel in- 
sights gained by such methodology, and we will continue 
to use this model system for this purpose here. Of course, 
this system is very interesting in its own right, too. Orig- 
inally, there was work on supported catalysts by Cant, 
Hicks and Lennon^, as well as on single crystals by Pe- 
den and Goodman 5 , showing that what they called Ru 
catalysts exhibited high activity for CO oxidation un- 
der reactant pressures of the order of atmospheres. Ex- 
tensive experimental and theoretical work has by now 
shown that in the corresponding reactive environments 
in fact Ru02 forms at the surfaceJiiSuLiSiiAii^, and ac- 
tuates the catalytic functioni&±ii 1 £i±&. Specifically, at 
the Ru(0001) model catalyst surface, it is an epitaxial 
RuO2(110) film that is formed, and although domain 
boundaries and steps are present, their influence on the 
catalytic function is not significant HJL Our theoretical 
modeling will therefore focus on the surface processes at 
a RuO2(H0) facet under steady-state conditions. The 
theoretical results will be compared with existing exper- 
imental data£*i&, and with results obtained previously 
within the computationally much less demanding "con- 
strained thermodynamics" approach 18 . 

Having already briefly communicated the highlights 
of this work 3 ., a large part of the present paper is de- 
voted to a detailed description of the employed method- 
ology. For this purpose we first establish in sections IIA 
and B the general theoretical framework for the com- 
bined DFT+kMC scheme, before addressing in section 
IIC the specific microscopic parameters calculated for the 
RuO2(110) system. The method is then employed in sec- 
tion III to analyze the surface composition and structure, 
showing that high catalytic activity in this system is inti- 
mately connected with a disordered and dynamic "phase" 
at the surface, in which the observable catalytic function 
results from the self-regulating action of several different 
elementary processes. 



II. THEORETICAL FRAMEWORK 

A. First-principles kinetic Monte Carlo simulations 

We aim at a first-principles statistical mechanics de- 
scription, which provides atomic-scale understanding of 
the steady-state catalytic function by explicitly consid- 
ering the detailed statistical interplay of all elementary 
processes, i.e. by fully accounting for correlations, fluc- 
tuations and spatial distributions. This requires consid- 
eration of a large enough system size, but due to the 
kinetic nature of steady-state catalysis also to follow the 
evolution over time scales that are long enough to ob- 



tain reliable averages. Particularly the latter is quite 
demanding, since elementary processes in heterogeneous 
catalysis are typically activated and thus rare. While in- 
dividual events occur on electronic to atomic time scales 
(1CP 15 — 10~ 12 sec), the time between consecutive events 
can be many orders of magnitude longer, possibly ren- 
dering it necessary to simulate up to seconds or more in 
order to arrive at meaningful conclusions about the effect 
of the statistical interplayi 

Kinetic Monte Carlo (kMC) algorithms suitably ad- 
dress this problem by providing a numerical solution 
to the Markovian master equation, which describes the 
dynamic system evolution by efficiently coarse-graining 
the molecular dynamics to the decisive rare events 
and properly averaging over the irrelevant short-time 
dynamics>i2i22i2ii22iS^ This allows to access time scales of 
the order of seconds or longer, even for mesoscopically- 
sized systems, while still retaining the full atomistic 
information^ A kMC simulation proceeds by generat- 
ing a sequence of system configurations. At each con- 
figuration, all possible elementary processes p and the 
rates r p at which they occur are evaluated. Appropri- 
ately weighted by these different rates one of the possi- 
ble processes is then executed randomly, leading to the 
next system configuration. In practice this is achieved by 
determining the total rate R = Ylp=i r p summing over 
all P possible processes, and executing process k, which 
fulfills the condition 

k fc-1 

5> p > Pl R > ^r p , (1) 
p— 1 p— l 

where p\ G [0, 1[ is a random number. This way, the 
kMC algorithm effectively simulates stochastic processes 
described by a Poisson distribution, and a direct and un- 
ambiguous relationship to real time is established by ad- 
vancing the clock by 



where pi S [0, 1[ is another random number. 23 

The crucial ingredients to a kMC simulation are there- 
fore the analysis and identification of all possibly rel- 
evant elementary processes at a given system configu- 
ration, and the determination of the associated rates. 
These rates are hitherto often either guessed or fitted 
such that the results of the kMC simulations match e.g. 
some existing experimental data. Modern first-principles 
kMC simulationa 25 i 26 i 27 i 2S i 29 use instead rates obtained 
by electronic structure theory calculations, ensuring that 
the parameters fed into the kMC simulation have a clear 
microscopic meaning. To keep the number of elemen- 
tary processes (and thus the number of required rates) 
tractable, it can be useful to suitably coarse-grain the 
system onto a lattice. This is the approach we will em- 
ploy here, and for the present case of CO oxidation at 
RuO2(H0) this then leads to an exhaustive list of 26 
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processes. We note, however, that for more complex sys- 
tems it may not be straightforward to identify a suitable 
lattice, or even if a lattice can be found, it may be diffi- 
cult to identify and consider all of the possible processes 
occurring on iti^ii 

B. Rates for elementary processes 

We now describe a consistent framework for the com- 
putations of the rates of elementary processes needed for 
a first-principles kMC simulation. Since the employed 
methodology is not restricted to the application to the 
RuO2(110) model catalyst, we describe the framework in 
this section in more general terms, and refer only occa- 
sionally to the implementation in the RuO2(110) system 
for illustration. 

Let us consider a discretized lattice model of a surface, 
where each surface unit-cell of area A uc is represented by 
st different site types (e.g. bridge and cus sites in the 
application to RuO2(H0) discussed below). This surface 
is exposed to a thermal gas composed of i different species 
with masses rrii. At the accuracy level relevant for our 
study, this gas is well described by ideal gas laws, and 
is then sufficiently characterized by the temperature T 
and partial pressures pi . In the context of heterogeneous 
catalysis, relevant process types occurring in the system 
are adsorption and desorption processes between the gas 
phase and the lattice, surface chemical reactions, as well 
as diffusion processes between lattice sites. 

1. Adsorption and desorption: General 

The adsorption rate of species i onto a free site of type 
st depends on the kinetic impingement onto the whole 
surface unit-cell, as well as a local sticking coefficient 
S s t.i{T), which governs which fraction of these impinging 
particles actually sticks to the given free site, 

r&(T, Pi ) = S st>i (T) PiA ™ , (3) 
V ZirrriikB-L 

where fcs is the Boltzmann constant. The impinging par- 
ticles from the thermal gas have initially a randomly dis- 
tributed lateral position of their center of gravity over 
the unit-cell, a random distribution over their possible 
internal degrees of freedom (e.g. molecular orientations) 
and Maxwell-Boltzmann distributed velocities. The lo- 
cal sticking coefficient represents thus a statistical aver- 
age over these degrees of freedom of the impinging gas 
phase particles of species i, and gives only the fraction 
that sticks to a given free site of the specified site type 
st. 

If there is only one site type per surface unit-cell, 
S s t,i{T) is identical to the often studied initial sticking 
coefficient, S 0} i(T)M^ The latter also describes a simi- 
lar statistical average of sticking probabilities, but does 
not distinguish between site types. It thus represents the 
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FIG. 1: Potential-energy surface (schematic) along the mini- 
mum energy path connecting the free gas phase and adsorbed 
state in surface site st. Shown is a case, where the adsorption 
process is exothermic (Ef, yB t,i < 0), but activated (AE^i > 0). 

average sticking to all sites of a given clean surface and 
we have 

So AT) = Y,S s t,i( T ) ■ ( 4 ) 

st 

If there are different adsorption site types in the surface 
unit-cell, knowledge of S ^(T) is not sufficient for a kMC 
simulation. For the latter we need to know specifically 
the rate with which particles adsorb in each of the differ- 
ent available site types, and the relevant quantity is then 
SstAT)- 

In the so-called hole model for adsorption^, one con- 
siders that only gas phase particles of species i with initial 
lateral position within a given area A st ^ around the site 
st have a chance to stick to the site. The local sticking 
coefficient will then contain a term ( ^ 3f '' ) < 1, reducing 
the overall impingement onto the whole surface unit-cell 
contained in eq. © to this active area. Note, that due to 
this term the overall adsorption rate r^^T,^) does also 
not depend on the actual choice of the surface unit-cell; 
only the gas phase particles impinging through the ac- 
tive area A stj i contribute to it. In this respect, the total 
sum of active areas can also not exceed the whole surface 
unit-cell area, since then more particles would be consid- 
ered than the actual overall impingement rate. Hence, 
the choice of active areas is additionally limited by the 
condition 

^ ] A s t,i 5-- A uc . (5) 

st 

In a classical picture, sticking occurs if the particle's 
initial velocity is high enough to surmount the barriers 
along the trajectory to the surface. If all particles of 
species i impinging within the active area A st ^ around 
the site type st had to pass the same maximum ener- 
getic barrier AE^^ > 0, the local sticking coefficient 
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would simply be given by ( A ^ t '' ) exp ^— A fc ^'' ^ ^ For 

a general form of the high-dimensional potential-energy 
surface (PES) seen by the impinging particles, this can 
be suitably generalized to 
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exp 



(6) 



Here, AE^f i denotes the maximum barrier height along 
the minimum energy path (MEP) connecting the gas 
phase and adsorbed state in the PES as schematically 
shown in Fig. ^ The factor /^(T) < 1 accounts for 
a further reduction in the sticking probability, if parti- 
cles with certain initial states (lateral positions over the 
unit-cell, internal degrees of freedom) are not efficiently 
steered along the MEP and are then possibly reflected by 
a higher barrier. 

Desorption of particles from occupied sites of the lat- 
tice represents the time-reversed process to adsorption, 
and as such has to fulfill detailed balance or microscopic 
reversibility. This principle ensures that the description 
of the free energy landscape for the forward and back- 
ward process between gas phase and adsorbed state is 
done consistently, and requires that the rate of adsorp- 
tion and desorption are related by 
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, (7) 



Here, AG s t^(T,pi) is the change in Gibbs free energy 
between the gas phase and the adsorbed state at site st. 
Neglecting the small piV term for the adsorbed stated, 
this is approximately given by the difference between the 
chemical potential fj, ga , St i(T,pi) of the particle in the gas 
phase and the free energy F sM (T) = - k B T\n{q v s 1 ^) 
in the adsorbed state. El°\ is hereby the total energy, 
and the partition function gj^ accounts for the vibra- 
tional degrees of freedom in the adsorbed state. 

The chemical potential of an ideal gas can be writ- 
ten as sum of total energy E^°l i and free energy con- 
tributions from the translational and internal modes 
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where g™* i is the partition function accounting for the 
internal degrees of freedom of the free gas phase parti- 
cle. With the chosen sign convention, [i gas ,i(T,pi) ap- 
proaches thus — oo in the limit of an infinitely dilute gas. 
For diatomic gas molecules like O2 and CO in the present 



case, AyU gaS) i(T,}3j) can be readily calculated from first- 
principles, yielding results that are, at all temperatures 
and pressures of interest to us here, virtually indistin- 
guishable from the experimental values listed in thermo- 
chemical tables^ 

Introducing the binding energy for the particle i at site 
st and with respect to its reference state in the gas phase, 
E b>st ,i = £*M-£gas,* ( with Eb,st,i < for exothermicity), 
cf. Fig. allows to rewrite eq. Q to 



,.ad 



(T,Pi 



„dcs 



(T) 



exp 



^ ACV.,,,i.7-.y; ; : 



— *c exp 
vib f 

Ist.i 



k B T 

A^gas,i(T,Pi) - E] 



b,st,i 



k R T 



(9) 



Equations ©, JSJl, and provide a consistent frame- 
work for the description of the kinetics of adsorption and 
corresponding desorption processes. The system-specific 
quantities required to fix the adsorption and desorption 
rates via these equations are f^^T), A st ^ and AE^ i 
describing the adsorption process, and Eb, s t,i and qj^ 
describing the bound state at the surface. The equations 
apply equally to processes that are unimolecular in ei- 
ther direction, or to processes that describe dissociative 
adsorption and in turn associative desorption. The pro- 
cesses can involve only one site of the surface lattice, as 
e.g. in the case of non-dissociative CO adsorption and 
desorption in the RuO2(110) example discussed below. 
Or they can involve more sites, possibly due to dissocia- 
tive adsorption or because a larger molecule blocks sev- 
eral sites upon adsorption. In the catalytic context, also 
surface chemical reactions involving direct desorption of 
the product may be described with the same framework 
of equations, viewing such processes then as dissocia- 
tive adsorption and associative desorption of the product 
molecule. 



2. Adsorption and desorption: Transition state theory 

The computationally most demanding quantity en- 
tering the hitherto established framework is the factor 
f%ti(T), accounting for a reduction in the sticking prob- 
ability due to those particles not efficiently steered along 
the MEP. An accurate calculation of /^(T) requires dy- 
namical simulations of a statistically meaningful number 
of trajectories of impinging particles, for which informa- 
tion about large parts of the underlying high-dimensional 
PES is required. 32 In the case of activated adsorption, 
an approximate f%fi(T) can also be determined through 
transition state theory (TST), relying on less PES in- 
formation. As described above, /f t d j(r) summarizes the 
dependence on all other degrees of freedom of the gas 
phase particle apart from the reaction coordinate along 
the MEP. In TST this is approximated by the fraction 
of correspondingly accessible states for thermalized par- 
ticles at the MEP barrier and in the initial gas phase, i.e. 
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Here, q^g^ st ^ is the partition function at the barrier 
(i.e. the transition state, TS), and we have separated 
the total partition function in the initial state conve- 
niently into contributions from the remaining two lateral 
translational degrees of freedom (?*asT :2D = ^uc 27Tm ^ BT 
over the surface unit-cell, and the internal degrees of free- 
dom <?g"g j of the free gas phase particle. The information 

about the PES required for an evaluation of /^ TST (T) 
is thus the location of the TS, as well as the local PES 
around it. 

Equation l|l(J|) introduces an approximation, since in 
contrast to the unspecified f^f i (T) the fraction of parti- 
tion functions accounts no longer for possible recrossings 
of the impinging particles after they have successfully sur- 
mounted the MEP barrier, and also limits the accessible 
states at the TS to those of an equilibrated particle. To- 
gether with the initially made assumption of a classical 
barrier, these are the three classic assumptions underly- 
ing TST—. Indeed, using / st | T T (T) in the adsorption 
rate in the detailed balance equation, eq. 10, yields an 
expression for the desorption rate familiar from TST 



T ) 



where 



Jst.i y 1 ) 



AE 



dcs 



exp 



k B T 



Jst,i K 1 I 



( Aihi „vib 



st.i) 



(11) 



(12) 



and we have realized that along the MEP, the adsorption 
barrier AE^^ and desorption barrier AE^\ are related 
through, cf. Fig. ^ 

E b , st ,i . (13) 



AEf™ = AE** 



3. Diffusion 

Diffusion processes of adsorbed particles on the lattice 
can also be treated within a TST framework. The rate 
for such hops from one site (of type st) to another site 
(of type st') is then given by 
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AEfl s st , i denotes the maximum barrier along the MEP 
between the two sites, <2 T g( st st i ^ the corresponding par- 
tition function at this barrier, and gj t lb the partition func- 
tion at the bound state in site st as before. 

Similar to the process types covered before, fulfillment 
of microscopic reversibility is also here a crucial issue. 
This is trivially given for diffusion between identical site 
types, since f^J^ 8 ' 1 '(T) and AEfl s st i are identical for 
the forward and backward jump. For diffusion between 
different site types, e.g. bridge to cus site and the re- 
versed process cus to bridge site in the RuO 2 (H0) ex- 
ample, this is no longer the case. If approximations are 
employed here, they have to be made consistently for 
both the forward and backward diffusion process. 



C. First-principles parameters for the CO 
oxidation over RuO2(110) 

1. Lattice model 

The rutile RuO2(110) surface exhibits a rectangular 
(lxl) unit-cell as shown in Fig. [3 In this orienta- 
tion, the bulk stacking sequence can be expressed as a 
series of 0-(Ru202)-0 trilayersii2i^4 Depending at which 
layer the sequence is truncated, three different (1 x 1) 
surface terminations result, and the middle panel of Fig. 
121 shows a topview of the surface atomic structure when 
the cut yields a top (RU2O2) layer. All oxygen atoms 
in this topmost (RU2O2) layer maintain their threefold 
bulk-like coordination to Ru atoms and are thus tightly 
bound. 18 Binding to the undercoordinated surface Ru 
atoms, on the other hand, gives rise to two prominent 
adsorption sites per surface unit-cell, namely a bridge 
(br) site bridging two surface Ru atoms and a so-called 
coordinatively unsaturated (cus) site atop of one sur- 
face Ru atom, cf. Fig. |3 The rutile bulk-stacking se- 
quence would be continued by oxygen atoms occupying 
first all bridge sites (leading to the so-called O br /- termi- 
nation, upper left panel in Fig. [2J and then all cus sites 
(O br /O cus termination, lower left panel), before another 
(Ru 2 2 ) plane would start a new sequence. The O br /- 
termination is routinely observed experimentally after 
high temperature anneals in ultra- high- vacuum (UHV), 
whereas the O br /O cus termination is most stable a high 
oxygen pressuresiii*i2i2i Extensive theoretical work over 
the last years has shown that O and CO adsorption at 
other sites of the surface is energetically significantly less 
favorable than adsorption at bridge and cus sitesAi^i 8 ^^, 
which is consistent with all presently available experi- 
mental data 7 i n i 37 i 38 . 39 i 40 i 41 i 42 . 

In view of these findings we choose the following model 
of the RuC>2(110) system as basis of our kMC simula- 
tions: The surface is described by a periodic lattice with 
a rectangular surface unit-cell, each cell containing one 
bridge and one cus site. The arrangement of these sites 
is such that along the [110] direction, i.e. the long side of 
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cus bridge 
cus site bridge site s j tes sites 




FIG. 2: Top view of the RuO2(110) surface showing the two prominent adsorption sites (bridge and cus), which continue 
the bulk-stacking sequence (central panel). When focusing on these two site types, the surface can be coarse-grained to the 
lattice model shown schematically to the right. Additionally shown are two perspective views to the left, exemplifying what 
the atomic structure of the surface looks like, if all bridge sites are occupied with oxygen atoms and the cus sites remain empty 
(O bl /- termination, top left panel), or if oxygen atoms occupy both site types (O br /O cus termination, bottom left panel). Ru 
= light, large spheres, O = dark, medium spheres. Atoms lying in deeper layers have been whitened in the topview for clarity. 



the rectangle, bridge and cus sites are alternating, while 
along the short direction of the rectangle (representing 
[001]) only like sites align. This leads to a characteristic 
stripe pattern with alternating rows of bridge and cus 
sites as schematized in the right panel of Fig. |2J Each of 
the sites in the lattice can be either empty, or be occu- 
pied by one O or CO. Since this comprises the possibility 
of bridge sites not being occupied with oxygen atoms, 
we note that our treatment includes thus implicitly the 
effect of the frequently discussed O surface vacancies in 
the UHV O br /- termination. 

In the context of catalytic CO oxidation, the kMC sim- 
ulations need to consider adsorption, desorption, diffu- 
sion and chemical reactions on this lattice. CO adsorp- 
tion into vacant cus or bridge sites is non-dissociative, 
while oxygen adsorption is dissociative and requires two 
vacant neighboring sites, i.e. a br-br, cus-cus, or br-cus 
pair, cf. Fig. [21 Together with the corresponding (time- 
reversed) desorption processes this leads already to 10 
different process types. Diffusion of adsorbed O or CO 
is modeled as hops to nearest neighbor lattice sites, and 
since adsorbatc diffusion can go br-to-br, br-to-cus, cus- 
to-cus, or cus-to-br, this adds another 8 processes. Prod- 
uct CO2 molecules may finally be formed by reaction of 
adsorbed O and CO in neighboring sites. There are there- 
fore four different reaction processes possible, namely 
O br +CO br , O br +CO cus , O cus +CO br and O cus +CO cus . 
At the temperatures of interest to our study the formed 
CO2 desorbs immediately^ , so that one reaction event 
is modeled as associative CO2 desorption and leads to 
two vacant sites on the lattice. In order to fulfill detailed 
balance, also the back reaction, i.e. dissociative CO2 ad- 
sorption, has then to be included, adding another four 
processes to the list. In total, our kMC simulations on 
this lattice model consider therefore 26 different elemen- 



tary processes, for each of which DFT calculations were 
carried out to provide the parameters needed to deter- 
mine the process rates via the equations derived in Sec- 
tion IIB. 



2. DFT setup 

The DFT calculations use the full-potential linear 
augmented plane wave approach (FP-LAPW)^^^ 
together with the generalized gradient approximation 
(GGA) for the exchange-correlation functional-?. The 
reference calculations for the free gas phase molecules, 
as well as all calculations addressing O2 adsorption and 
desorption barriers were done in a spin-polarized mode. 
Extensive test calculations confirmed that in all adsorbed 
configurations, including surface diffusion, the spin is 
quenched at the metallic RuO2(110) surface, and the cor- 
responding computations were carried out in a non-spin- 
polarized way. The RuO2(110) surface is modeled in a 
supercell geometry, employing a symmetric slab consist- 
ing of three rutile 0-(Ru202)-0 trilayers. All atomic 
positions within the outermost trilayer were fully relaxed 
and a vacuum region of w 11 A ensures the decoupling 
of the surfaces of consecutive slabs. 

The FP-LAPW basis set parameters are: -R^t =1-8 
bohr, i?MT =1-1 bohr, -R£t =1-0 bohr, wave function 
expansion inside the muffin tins up to Z™a X = 12, poten- 
tial expansion up to lf°^ — 4. For the RuO2(110) slabs 
the (1 x 1) Brillouin zone (BZ) integration was performed 
using a (5 x 10 x 1) Monkhorst-Pack grid with 50 (15) 
k-points in the full (irreducible) part of the BZ. To ob- 
tain the same sampling of the reciprocal space for bigger 
surface cells, this number is reduced accordingly. The 
energy cutoff for the plane wave representation in the in- 
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terstitial region between the muffin tin spheres was i?™ f ax 
= 20 Ry for the wave functions and = 169 Ry for 

the potential. 

The computational setup is thus exactly the same 
as the one employed in our previous work on 
RuO 2 (110)i 18 i 34 Its numerical accuracy has already been 
detailed there, and with respect to the present study it 
is characterized by a ±150 meV uncertainty in the bind- 
ing energies and activation barriers, and a ±50 meV un- 
certainty in relative binding energy differences for the 
same species (e.g. at different coverages). Although 
already this setup led to a massive computational de- 
mand, the absolute numbers are therefore not well con- 
verged. We confirmed, however, that the correct ener- 
getic ordering is obtained 48 , which is what largely deter- 
mines the results discussed in this work. The effect of 
the uncertainty in the kMC rates due to the errors in 
the energetics will be further discussed in section IIIC. 
This then also comprises the additional uncertainty intro- 
duced by the more basic deficiency of density-functional 
theory, namely the approximate nature of the employed 
exchange-correlation functional, here the GGA. 



3. Adsorption, Desorption and Reaction 

Within this setup we first address processes, which 
change the number of oxygen atoms or CO molecules ad- 
sorbed at lattice sites. This naturally includes adsorption 
and desorption of the reactants, but as mentioned above 
also reaction processes leading to CO2 formation, which 
are treated as associative CO2 desorption. The energet- 
ics along the MEP summarized in Fig. 2] play a decisive 
role in determining the rates of these processes. In gen- 
eral, all corresponding quantities appearing in eq. i|13|) 
depend not only on the site types actually involved in the 
process, but also on the local environment around them, 
i.e. on the occupation of nearby lattice sites. Formally 
such adsorbate-adsorbate interactions can be taken into 
account through a lattice gas Hamiltonian (LGH), ex- 
panding the total energy in terms of on-site energies and 
lateral interactions of pair or higher many-body typei 

In order to assess the variations in the process energet- 
ics that arise due to lateral interactions we concentrate 
first on the binding energy, i.e. the minimum of the MEP. 
In the sense of a LGH, lateral interactions affecting this 
quantity may be classified into interactions between like 
(i.e. br-br, cus-cus) and unlike (i.e. br-cus) sites, as well 
as into interactions between like (i.e. O-O, CO-CO) and 
unlike (i.e. O-CO) species. Starting with the interac- 
tions between like sites and like species, we computed a 
series of surface configurations involving always only one 
surface species and occupation of only one site type at 
various concentrations. Specifically, we employed (lxl), 
(1 x 2) and (1 x 3) surface unit-cells to describe atomic 
arrangements, in which every, every second or only every 
third site along the rows in the [001] direction, cf. Fig. |3 
is occupied. Calculations in (2 x 1) cells allowed accord- 



TABLE I: Binding energies, Eh (in eV), in (1 x 1) surface 
configurations involving only one species (O or CO) and occu- 
pation of all sites of one site type (br or cus), i.e. a coverage of 
1 monolayer (ML) . Using surface unit-cells of varying size, this 
total coverage is varied between 1/3 and 1 monolayer (ML), 
and the resulting variations in binding energy with respect to 
the value at 1 ML are given. 
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across rows 
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(1x2) 


(1x3) 
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\ ML 


i ML 


I ML 


Qbr 


-2.44 


+0.11 


+0.09 


+0.01 


QCUS 


-1.08 


-0.06 


-0.04 


-0.06 


co br 


-1.64 


-0.02 


-0.02 


+0.08 


co cus 


-1.31 


-0.07 


-0.07 


+0.07 



ingly to change the concentration in the [110] direction 
across the rows. The resulting binding energies are sum- 
marized in Table [I] They show only minor variations, 
indicating overall rather weak lateral interactions of O br - 
O br , O cus -O cus , CO br -CO br , and CO cus -CO cus type. 

We interpret this finding as a consequence of the rather 
open crystal structure of the RuO2(110) surface. For a 
bridge site, already a nearest neighbor bridge site across 
the rows and a second nearest neighbor bridge site along 
the rows is more than 6 A away, cf. Fig. [21 and neither 
of these sites involves a coordination to surface atoms al- 
ready involved in the direct bonding to the given bridge 
site. For a cus site and its neighboring cus sites, this sit- 
uation is exactly the same. A sharing of surface atoms at 
close distance occurs only upon occupation of the near- 
est neighbor bridge site along the bridge rows, explaining 
the slightly larger variation of the binding energy for a 
surface atomic arrangement in which oxygen atoms oc- 
cupy such neighboring sites, cf. Tabic [I] It is interest- 
ing to note that the corresponding binding energy varia- 
tion is masked for the case of adsorbed CO br molecules, 
since there the bonding geometry changes between a sym- 
metric and an asymmetric configuration, depending on 
whether neighboring bridge sites are occupied or not. 
This has been discussed extensively in Refs. andl42l 
and our DFT results concerning the structural and ener- 
getic details are in full agreement with the data reported 
there. 

The variations in the binding energies due to lateral 
interactions between like species and like sites are always 
well within ±150 meV compared to the value in a (1 x 1) 
arrangement. Similar findings were obtained in calcula- 
tions for atomic arrangements involving again only occu- 
pation of one site type (either bridge or cus) , but having 
O and CO simultaneously present at the surface. Again, 
the binding energy variations were well within ±150 meV, 
indicating that also interactions of the type O br -CO br 
and O cus -CO cus are rather small. This leaves only inter- 
actions between bridge and cus sites, both between like 
and unlike species. Calculations with sparse atomic ar- 
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TABLE II: Binding energies, Eb (in eV), of O and CO surface 
species adsorbed in a (1 x 1) configuration and with varying 
occupations at the other site type within the surface unit cell. 





occupation at nearest neighbor bridge site 
empty O br CO br 


QCUS 

co cus 


-1.08 -0.99 -1.01 
-1.31 -1.26 -1.26 




occupation at nearest neighbor cus site 
empty O cus CO cus 


Qbr 

co br 


-2.44 -2.23 -2.37 
-1.64 -1.57 -1.58 



rangements in (1 x 2) and (1 x 3) surface unit-cells showed 
in this case that br-cus interactions between second and 
third nearest neighbor sites at distances larger than 4.4 A 
are negligibly small. The binding energy variations be- 
come only somewhat more appreciable, when nearest 
neighboring bridge and cus sites are involved. This is 
illustrated in Table ITT1 for arrangements with (lxl) pe- 
riodicity, showing that primarily adsorbed oxygen atoms 
exhibit some variation in the bond strength, depending 
on whether the other site type within the surface unit-cell 
is either empty or occupied by O or CO. 

In summary, our calculations reveal that lateral in- 
teractions, affecting the binding of O and CO to the 
RuO2(110) surface, are small. The only significant in- 
teractions are with nearest neighbor lattice sites at an 
around 3 A distance, cf. Fig. [21 Along the rows in 
the [001] direction this is a like site (i.e. br-br or cus- 
cus), whereas across the rows in the [110] direction this 
is the interaction with the next unlike site (i.e. br-cus). 
We expect such interactions to play a role in stabiliz- 
ing ordered adsorbate superstructures at the surface at 
low temperatures. With E b (O cus ) = -1.0 ± 0.15 eV, 
E b (O hl ) = -2.3 ± 0.15 eV, £ b (CO cus ) = -1.3 ± 0.15 eV, 
and E b (CO hr ) = -1.6 ± 0.15 eV in all calculated test ar- 
rangements, we recognize, however, that the variations 
are still rather small compared to the other energies in 
this system. For the present study, we therefore ne- 
glect adsorbate-adsorbate interactions and employ the 
just listed four E b values to describe the binding of O 
and CO at the two sites throughout. 

Considering the results obtained for the binding ener- 
gies, we do not expect the other parts of the MEP to 
be much more affected by lateral interactions. All cal- 
culations of energy barriers along adsorption, desorption 
and reaction pathways were correspondingly carried out 
using (1x2) surface unit-cells, this way ensuring that 
at the relevant transition state geometries the distances 
to all periodic images are always larger than 6 A. Start- 
ing the description of such barrier calculations with the 
adsorption and desorption of the reactants, we first ad- 
dress the associative oxygen desorption from two neigh- 
boring cus sites. For this we computed the PES, when 
lifting the associating molecule from the O br /O cus termi- 
nation, i.e. using the equal height of both O cus atoms 




2 bondlength (A) 



FIG. 3: Potential-energy surface (PES) for associative oxygen 
desorption from two neighboring cus sites into an O2 molecule 
with orientation parallel to the surface and aligned along the 
[001] trench (as schematically shown in the inset). The height 
above the topmost (RU2O2) surface and the distance between 
the two O atoms are used as coordinates of the plot. For 
every (height, bondlength) of the parallel O2 molecule, all 
atomic positions of the oxide substrate are fully relaxed. Ac- 
tually calculated points are indicated, and the energy zero 
corresponds to a free O2 molecule far above the surface. 

and the distance betweem them as reaction coordinates, 
while minimizing the energy at each point with respect 
to all remaining degrees of freedom. As apparent from 
Fig. |3 there is a shallow molecular chemisorption well 
at a height of 1.8 A, which has already been described in 
previous DFT calculations^!. Above this well the energy 
rises, however, smoothly to the level of a free gas phase 
molecule, reflecting that for this molecular orientation 
dissociation is non- activated. 

Equivalent findings were also obtained when lifting a 
pair of O atoms from neighboring bridge-bridge sites and 
from neighboring bridge-cus sites, as well as when lift- 
ing up an upright CO molecule from a bridge and a cus 
site, i.e. none of these cases showed an activation barrier. 
While these calculations for fixed molecular orientation 
and lateral position over the unit-cell do, of course, not 
account for the variations in the full high-dimensional 
PES for adsorption and desorption, they still show that 
the maximum barrier height along the MEP connecting 
the gas phase and adsorbed state must be zero. As dis- 
cussed in section IIB.l we correspondingly use 

A pad a ™d a pad 

^^O.cus-cus — Lya O, br-br — ^-^O.cus-br — 

A££ d 0jCUS = AE& M = (16) 

in the determination of the local sticking coefficient 
through eq. ©. This leaves then only the quantities 
fsti{T) and A st ^ to entirely determine the rates of all 
5 possible adsorption processes of reactants, which in 
turn also determines the corresponding desorption rates 
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through eq. ©, using gj* = 1. 

The quantity /^(T) accounts for a reduction of the 
sticking probability, if particles with certain initial states 
(lateral positions over the unit-cell, internal degrees of 
freedom) are not steered along the MEP and are then po- 
tentially reflected by a higher barrier. In the temperature 
range of interest to the present study (~ 300 — 600 K), 
impinging thermal molecules are still comparably slow, 
and we expect them to be efficiently steered along non- 
activated pathways in the high-dimensional PES. As de- 
scribed above, an accurate determination of the there- 
fore expected value of /^(T) m 1 would require com- 
putationally highly demanding dynamical simulations of 
a statistically relevant number of trajectories for the im- 
pinging molecules^ As we will discuss in section IIIC, 
uncertainties in the individual rates by a factor of 10 to 
100 do not affect the reported conclusions. Since f%ti(T) 
enters the rates only linearly, we therefore simply set it 
to unity, and with the same argument, partition the ac- 
tive areas for adsorption equally over the various sites. 
Since there are four ways for an O2 to dissociate into site 
pairs in a (1 x 2) unit-cell, this leads to Ao, C us-cus = 
^o,br-br = ^o,cus-br = 1/2A UC , and with the same ar- 
gument we obtain ^4co,cus = ^4co,br = 1/2 A uc for CO 
adsorption. With this choice of prefactors the sum of lo- 
cal sticking coefficients over all site types in the surface 
unit-cell has the correct upper bound, and through the 
use of eq. (JSJ to determine the corresponding desorption 
rates, detailed balance is also ensured. 

Having determined the rates of the 10 possible adsorp- 
tion and desorption processes of the reactants, we are 
left with surface reactions as the final event type chang- 
ing the number of O or CO adsorbed at lattice sites. 
Similar to the procedure employed for the desorption 
events we evaluate the energetics by mapping out the 
PES along suitable reaction coordinates, fully relaxing 
all remaining degrees of freedom. For two of the four 
possible reaction processes, this was already described 
in detail in a preceding publication^, and exactly the 
same methodology was now used to locate the remain- 
ing two transition states. The resulting reaction barri- 
ers (or equivalently associative CO2 desorption barriers) 
are: AE^ b s I+CQbI = 1.54 eV, AEg£ +GOou . = 1.25 eV, 
A^L+cqcu, = 0.89 eV, and AE^ +cobt = 0.76 eV. 
Since in this case all four reaction processes have clear de- 
fined transition states, it is useful to resort to eq. (fTTJl to 
fix the process rates. For this, we employ ff^l' TST = 0.5, 
which includes a value of A s t.i — 1/2 A uc , again simply 
equally partitioning the active areas for the correspond- 
ing time-reversed dissociative CO2 adsorption processes. 
Effectively, we therefore set g^s (st ») = Qst?i > CI - ec l- > 
and a more rigorous approach would be to explicitly 
evaluate the computationally more demanding partition 
function q^s^st i) a ^ the barrier. Test calculations in the 
harmonic approximation indicate that the uncertainty in- 
troduced by the employed approximation translates into 
variations of the desorption rate by a factor of 10 (at most 



TABLE III: Binding energies, Eh, for CO and O (with respect 
to (1/2)02) at bridge and cus sites, cf. Fig. |21 and diffusion 
energy barriers, A_E dlff , to neighboring bridge and cus sites, 
as used in the kMC simulations. The desorption barriers are 
given for unimolecular and for associative desorption with ei- 
ther O cus or O br . This includes therefore surface reactions 
forming CO2, which are considered as associative desorption 
of an adsorbed O and CO pair. All values are in eV. 
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100), i.e. an uncertainty that is of the same order as the 
one in the determination of the reactant adsorption and 
desorption process rates. 

Together with the reactant binding energies, the des- 
orption barriers can be used to obtain the activation bar- 
riers for the corresponding dissociative CO2 adsorption 
processes. This yields in all four cases sizable barriers, 
which when input into eq. I© lead to vanishingly small 
values below 10 -12 for the local sticking coefficient at all 
temperatures of interest. Regardless of the exact CO2 
partial pressures building up over the operating catalyst 
surface, dissociative CO2 readsorption into vacant sites 
at the surface will therefore be completely negligible com- 
pared to the non-activated adsorption of oxygen and CO. 



4- Diffusion 

With 18 process rates determined, we now turn to the 
8 diffusion processes. In line with the computations de- 
tailed above, the corresponding barriers were obtained 
using (1x2) unit-cells and mapping the PES along the 
high-symmetry line connecting the initial and final site 
of the diffusion process, while fully relaxing all remain- 
ing degrees of freedom. Test calculations checking on 
the effect of lateral interactions with the periodic images 
were performed in suitably extended surface unit-cells, 
namely (1x3) unit-cells for diffusion processes along the 
trenches (i.e. cus-to-cus and br-to-br), and (2 x 2) for 
diffusion across the trenches (i.e. cus-to-br and br-to- 
cus). The observed < ±80meV variations with respect 
to the barriers obtained in the standard (1x2) cells 
reflect just as in the case of the binding energies that 
lateral interactions with atoms beyond 6 A distance are 
insignificant at this surface. Since the binding energy 
calculations indicated slightly larger lateral interactions 
between nearest-neighbor bridge and cus sites, all diffu- 
sion processes along the trenches were furthermore cal- 
culated with different occupations of the corresponding 
other sites, e.g. the oxygen cus-to-cus barrier was com- 
puted with O br or CO br present at the neighboring bridge 
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sites, or with the bridge sites empty. Again, the bar- 
rier variations were with ±130meV only modest. There- 
fore neglecting adsorbate-adsorbate interactions, we em- 
ployed the energetic barriers listed in Table [ffi] through- 
out to describe the eight possible types of diffusion events 
on the lattice. For completeness Table Hill also compiles 
the energetic barriers for all other processes as used in 
the kMC simulations. 

To fully fix the diffusion process rates in eq. (|14|) . we 
use fsf S s '^ ST = 1, i.e. like in the surface reaction case we 
again approximate 9xs(st st< i) ~ Qst?i- From test calcu- 
lations in the harmonic TST approximation we estimate 
that this has an uncertainty of a factor of 10. We will see 
below that the ensuing (already quite low) uncertainty 
in the corresponding rates is for this process class par- 
ticularly inconsequential. Due to the high barriers, most 
surface diffusion processes have such low rates that they 
occur only very rarely anyway, and under the catalyti- 
cally most relevant high-pressure conditions virtually all 
surface sites are occupied, thereby additionally blocking 
the remaining diffusion events. 

D. kMC setup for steady-state catalysis 

Using the now completely determined 26 process 
rates, kinetic Monte Carlo runs were performed for 
(T, PgOtPo 2 ) conditions covering the temperature range 
from T = 300 K to T = 800 K and partial pressures from 
10~ 15 to 10 5 atmospheres. In each simulation for fixed 
(T, PcO)Po 2 )j the evolution of the system from an arbi- 
trary initial lattice configuration is first calculated un- 
til steady-state conditions are reached, and the targeted 
average surface quantities are thereafter evaluated over 
sufficiently long time spans to reach a 1 % statistical pre- 
cision. The simulations were routinely repeated using 
other initial lattice configurations to verify that the true 
dynamic steady-state was reached. Besides the average 
surface populations, this kind of statistical analysis yields 
e.g. the average frequencies of the various elementary 
processes, and through the average frequency of all sur- 
face reaction events then also the absolute TOF under 
the chosen environmental conditions. 

The procedure is illustrated in Fig. 0] where the ac- 
tual and average surface populations of O br , O cus , CO br 
and CO cus are plotted, when starting from a fully oxygen 
covered surface (O br /O cus termination). Despite notable 
fluctuations, which we will discuss in more detail below, 
well defined average values for all surface species are fi- 
nally obtained. We note that although the individual el- 
ementary process dynamics takes place on a picosecond 
time scale, the "induction period" until the steady-state 
populations are reached is of the order of a tenth of a sec- 
ond. At the lower temperatures of interest to our study, 
namely 350 K, this becomes even more pronounced. Due 
to the decelerated process rates, the real times spanned in 
the simulations (at an equivalent number of kMC steps) 
can then be orders of magnitude longer. Only the efficient 
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FIG. 4: (Color online) Time evolution of the site occupation 
by O and CO of the two prominent adsorption sites, bridge 
and cus, cf. Fig. [5] The temperature and pressure conditions 
chosen (T = 600 K, pco = 20 atm, po 2 — 1 atm) correspond 
to optimum catalytic performance. Under these conditions 
kinetics builds up a steady-state surface population in which 
O and CO compete for either site type at the surface, as 
reflected by the strong fluctuations in the site occupations. 
Note the extended time scale, also for the "induction period" , 
until the steady-state populations are reached when starting 
from a purely oxygen covered surface. 

time coarse-graining underlying kMC algorithms makes 
it possible to reach such time scales, while still account- 
ing for the full atomic-scale correlations, fluctuations and 
spatial distributions. 

The simulations were carried out on a system with 
(20 x 20) surface sites (200 bridge and 200 cus sites) and 
periodic boundary conditions. Test runs for a (30 x 30) 
system gave identical results, and even then the compu- 
tational effort of the kMC simulations themselves was 
still insignificant compared to the first-principles calcu- 
lations required to evaluate the process rates. Before 
commencing with the real simulations, a useful check of 
the detailed balance of our setup was furthermore to run 
simulations with all surface reaction events switched off. 
The DFT-kMC approach is then equivalent to the com- 
putationally much less demanding "constrained equilib- 
rium" approach of atomistic thermodynamicsi^ and we 
indeed reproduce exactly our earlier results, which were 
based on a sub-set of the present DFT data—^i. The re- 
sulting average surface populations are shown in Fig. [S] 
to facilitate direct comparison with the kinetic data and 
are discussed in more detail below. 



III. RESULTS 

A. Steady-state surface composition and 
stoichiometry 

We start by using our first-principles kMC simula- 
tions to obtain the steady-state average surface popu- 



11 



p @ 600K (atm) 

10 15 10 10 10 5 1 10 +s 




p @ 350K (atm) 

FIG. 5: (Color online) Steady-state surface structures ob- 
tained by first-principles kMC simulations, where the surface 
reaction events are switched off (see text). In all non- white 
areas, the average site occupation is dominated (> 90%) by 
one species, i.e. either O, CO or empty sites ( — ). The la- 
bels correspondingly characterize the surface populations by 
indicating this majority species at the bridge and cus sites. 
The unlabeled small region just appearing in the upper left 
edge of the figure corresponds to a fully CO-covered surface, 
CO br /CO cus . Note, that in the thermodynamic case the sta- 
bility regions of the various surface phases depend only on 
the gas phase chemical potentials (i — 0,CO), leading to a 
simple scaling of the pressure axes at different temperatures, 
cf. eq. ©. The only exception to this is the width of the 
white coexistence regions, which increases with temperature. 
The one shown here corresponds to 600 K. Above the dashed 
line bulk Ru02 is thermodynamically unstable against CO- 
induced decomposition (see text). 



lations in a wide range of (T,pco,Po 2 ) conditions that 
are of interest in the catalytic context. For this sys- 
tem this is specifically the temperature window between 
~ 350 — 600 K, where high catalytic activity has been re- 
ported experimentally for varying pressures ranging from 
UHV up to technologically relevant values of the order 
of atmospheres^ii*i&. Figure summarizes the results 
for T = 600 K and varying O2 and CO partial pres- 
sures. Under most pressure conditions, the site occu- 
pation is found to be dominated by more than 90 % by 
one species (O or CO) or vacant sites (— ). This gives 
rise to four defined steady-state surface structures. For 
the lowest O2 and CO partial pressures this is the stoi- 
chiometric O br /— termination, where virtually all bridge 
sites are occupied by oxygen atoms and all cus sites are 
empty. Higher O2 partial pressures stabilize the O-rich 
O br /O cus termination, where now also all cus sites are 
occupied by O atoms. Both terminations are routinely 
observed experimentally under UHV conditions, the prior 
directly after high-temperature annealing and the latter 
after oxygen post-exposureiii*21 At increased CO par- 
tial pressures, we find two additional kinds of structures, 
CO br /- and CO b 7CO cus , both of which have recently 
also been characterized in UHV following different prepa- 



ration recipesaJaaiiaa. 

The transition from one structure to the other with 
changing pressures is not abrupt, but occurs over a finite 
pressure range. Under such conditions, more than one 
species (and/or vacancies) occupies a significant fraction 
(> 10%) of either surface site type, and we delimit such 
coexistence regions in Fig. El by dotted lines. If the sys- 
tem were in equilibrium with the surrounding gas phase, 
such coexistence regions between stable surface phases 
could be described by evaluating the configurational en- 
tropy, and in the present case due to the absence of lat- 
eral interactions this would simply be Langmuir adsorp- 
tion isotherms 49 . In fact, when switching off the reaction 
events in the kMC simulations the simulated average oc- 
cupations in coexistence regions are indeed identical to 
those resulting from the appropriate analytical Langmuir 
expression, since the kMC algorithm automatically ac- 
counts for the configurational entropy. Even with the 
reactions switched on, Langmuirian average occupations 
are still found in the coexistence regions between O br /— 
and O b 70 cus , as well as between O br /- and CO br /- 
at the lowest partial pressures. This situation changes, 
however, for the higher partial pressures, and particu- 
larly in the coexistence region between O br /O cus and 
CO br /CO cus the kinetics of the on-going surface chemical 
reactions leads to average site occupations, which can no 
longer be described by simple isotherms, and where the 
actual occupations exhibit complex spatial distributions 
and undergo large fluctuations. Since these conditions 
correspond also to the catalytically most active state of 
the surface, we will further analyze them in more detail 
in the section IIIB below. 

Before studying the full kinetics it is useful to com- 
pare the steady-state populations with those resulting 
from our previous ab initio atomistic thermodynamics 
investigations, employing the "constrained equilibrium" 
concept l,1 §. In these calculations we neglect the effect 
of the on-going catalytic reactions at the surface and as 
described in section IID above the resulting average sur- 
face populations (now identically computed by the kMC 
approach with the reaction processes switched off) are 
shown in Fig. When the catalytic reactions are 
switched on, we obtain the results shown in Fig. 
Comparing the two plots, similarities, but also the ex- 
pected noticeable differences under some environmental 
conditions can be discerned. The differences affect most 
prominently the presence of oxygen atoms at the bridge 
sites, where they are much more strongly bound than 
CO. For the thermodynamic approach only the ratio of 
adsorption to desorption matters, and due to the ensuing 
very low desorption rate, O br is correspondingly stabi- 
lized even when there is much more CO in the gas-phase 
than O2 (left upper part of Fig. |5J. The surface reac- 
tions, on the other hand, are very efficient for removing 
this O br species. As result, under most CO-rich con- 
ditions in the gas-phase oxygen is faster consumed by 
the reactions than it can be replenished from the gas- 
phase. The kMC simulations, as they account for this 
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FIG. 6: (Color online) Left panel: Steady-state surface structures obtained by first-principles kMC simulations at T = 600 K 
(see Fig. Qjjfor an explanation of the labeling). Right panel: Map of the corresponding turn-over frequencies (TOFs) in site _1 s _1 
(10 15 cm~ 2 s _1 ): White areas have a TOF < 10 -4 site _1 s _1 (< 10 11 cm'V 1 ), and each increasing gray level represents one or- 
der of magnitude higher activity. Thus, the darkest gray region corresponds to TOFs above 10 2 mol. site _1 s _1 (10 17 cm'V 1 ), 
while in a narrow (pco,Po 2 ) region the TOFs actually peak over 10 3 site^s" 1 (10 18 cm s ). The plots are based on ~ 400 
kMC simulations for different (pco,f>o 2 ) conditions. 



effect, then yield a much lower average surface concen- 
tration of O br than the thermodynamic treatment, and as 
a consequence show an extended stability range of sur- 
face structures with CO br at the surface (CO br /— and 
CO br /CO cus regions). It is particularly interesting to 
notice, that the resulting CO br /— region, consisting of 
only adsorbed CO at bridge sites, does not exist in the 
thermodynamic phase diagram at all, cf. Fig. [S] The cor- 
responding CO br /— "phase" (hatched region in Fig. |BJ) is 
thus a surface termination with defined atomic structure 
and composition that is entirely stabilized by the kinetics 
of the thermodynamically open catalytic system. 

The stability of this CO br /- and the CO br /CO cus 
structure in the most CO-rich environments has to be 
considered with caution though. In such environments 
one would expect a reduction of the whole oxide to pure 
Ru metal. In our kMC simulations, the two CO-covered 
surface terminations persist instead, because the simula- 
tions are restricted to the kinetics involving bridge and 
cus surface sites. This restriction is justified under most 
steady-state conditions, since the creation of vacancies in 
deeper layers, initiated by reacting off the tightly-bound 
threefold coordinated lattice oxygen atoms, is extremely 
rare. Even if such a vacancy would be created, it would 
be rapidly refilled by surface oxygen atoms. However, 
under the most CO-rich conditions and the resulting CO- 
covered surface terminations (upper left part of Fig. 0) 
such a rapid refilling will no longer occur and the ad- 
sorbed CO molecules will eventually reduce the R11O2 
catalyst to pure Ru metal. Under which environmen- 
tal conditions this instability of the bulk oxide against 
CO-induced decomposition occurs, can roughly be esti- 
mated through thermodynamic considerations^ and is 
indicated by the dashed line in Fig. [SJ As can be read 
off from the pressure scales in Fig. [SJ catalytically rele- 
vant environments are, however, quite far away from this 



instability limit. Although the steady-state kinetics re- 
duces the O concentration at the surface compared to 
this "constrained equilibrium" description, we still find a 
significant presence of surface oxygen atoms, cf. Fig. [5] 
which can rapidly anneal eventually created lattice oxy- 
gen vacancies. In all catalytically relevant conditions, 
only a small concentration of lattice oxygen vacancies 
will correspondingly be present during steady-state op- 
eration, and it is the kinetics of the bridge and cus site 
atoms and molecules which will entirely determine the 
catalytic function. 

Another important aspect that differs from the "con- 
strained thermodynamic equilibrium" results is a more 
complex dependence on temperature. In the thermody- 
namic approach the only ruling quantities are the gas 
phase chemical potentials, which means that for varying 
temperatures equivalent surface conditions are obtained 
at partial pressures which scale according to eq. (JSJl • We 
can therefore e.g. summarize the stability regions of the 
various surface structures at 350 K and at 600 K in Fig. [S] 
in the same range of chemical potentials by simply indi- 
cating two different pressure scales. The only notable ex- 
ception arises due to the varying configurational entropy, 
which leads to an increasing width of the white coexis- 
tence regions with temperature. When accounting for 
the kinetics, the situation becomes more complex. This 
is exemplified by Fig. which displays the average sur- 
face populations obtained in kMC simulations at 350 K 
and varying O2 and CO partial pressures. The pressure 
range shown in Fig. [7| corresponds to the same range of 
chemical potentials as the one in Fig. but due to the 
kinetics of the on-going surface reactions some differences 
between the results at 600 K and 350 K can now be dis- 
cerned. Scaling is still largely present though, since the 
change from oxygen- to CO-covered surface still occurs 
at roughly the same chemical potentials, indicating that 
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FIG. 7: (Color online) Left panel: Steady-state surface structures obtained by first-principles kMC simulations at T = 350 K 
(see Fig. |S]for an explanation of the labeling). Right panel: Map of the corresponding turn-over frequencies (TOFs) in site _1 s _1 
(10 15 cm~ 2 s _1 ): White areas have a TOF < 10 -12 site _1 s _1 (< 10 3 cm _2 s _1 ), and each increasing gray level represents one 
order of magnitude higher activity. Thus, the black region corresponds to TOFs above 10~ 5 site _1 s _1 (10 10 cm _2 s~ 1 ), while 
in a narrow (pco,Po 2 ) region the TOFs actually peak over 10 -3 site _1 s _1 (10 12 cm _2 s _1 ). The plots are based on « 400 
kMC simulations for different (pco,Po 2 ) conditions. 



the steady-state situation is mostly quite close to a "con- 
strained equilibrium" 18 . At the lowest partial pressures, 
the transitions are again Langmuir-like, although now of 
course with a reduced width corresponding to the lower 
temperature. Most important is that the surface compo- 
sition with the largest fluctuations is again found for the 
white coexistence region between the stability region of 
the O br /O cus and CO br /CO cus terminations, represent- 
ing also at this lower temperature the catalytically most 
active state. 



B. Catalytically active conditions 

To quantify the catalytic activity under the various 
conditions, the absolute TOFs were evaluated in the kMC 
simulations by analyzing the average occurrence of all 
surface reaction events. The results in the same range of 
partial pressures as discussed before are shown in Figs. 
El and [7] for 600 K and 350 K, respectively. As already 
mentioned above, the catalytic activity is at both tem- 
peratures narrowly peaked around the coexistence line 
between the O br /O cus and CO br /CO cus terminations, a 
region which we will henceforth shortly denote by the 
"active region" of the catalyst. Due to the afore discussed 
pressure scaling, this highest activity occurs at 600 K at 
technologically more relevant pressures of the order of at- 
mospheres, and at 350 K accordingly at pressures typical 
for UHV conditions. At 350 K the TOFs are significantly 
smaller than at 600 K. For both temperatures, the calcu- 
lated TOFs are in quantitative agreement with existing 
experimental data£*i&, both with respect to the absolute 
magnitude and their (T, po 2 iPco) dependence. 

A maybe even more remarkable aspect is revealed, 
when analyzing the contribution of the four possible sur- 



face reaction processes to the overall reactivity. As sum- 
marized in Tablc lITTl thcse four processes have quite differ- 
ent energy barriers, leading to large differences in the in- 
dividual rates. It is interesting to note that the energetic 
order among these four reactions may be rationalized on 
the basis of the initial bound states, since higher barriers 
are observed, when more strongly bound adsorbates are 
involved (in particular the strongly bound O bl species). 
On the basis of the energetic data, the lowest barrier 
process O cus + CO br involving the weaker bound surface 
species would have correspondingly appeared most rele- 
vant for the catalysis. Yet, in the narrowly peaked region 
of highest activity, its contribution is in fact only little. 
For instance, at po 2 = 1 atm and pco = 20 atm, i.e. 
within the most "active region" at T = 600 K, cf. Fig. El 
only 23% of the total computed TOF of 4.6 10 3 site -1 s -1 
(4.6 10 18 cm -2 s- 1 ) are due to this O cus + CO br reac- 
tion. With 76% it is instead the O cus + CO cus reaction, 
which dominates the overall activity, while a remaining 
1% is due to the O br + CO cus reaction. The highest 
barrier reation O br + CO br — > CO2 is insignificant. 

Although the quantitative shares of the different re- 
action processes exhibit some variations, this situation 
holds qualitatively for all environmental conditions corre- 
sponding to the "active region" at both 600 K and 350 K: 
The overall TOFs are in all cases dominated by the O cus 
+ CO cus reaction, although it exhibits a less favorable 
energy barrier compared to the O cus + CO br reaction. 
This somewhat surprising finding becomes comprehen- 
sible, when analyzing the steady-state surface popula- 
tions in more detail. Fig. shows a representative 
snapshot at 600 K, for peak activity pressure conditions 
(Po 2 — 1 atm, pco — 20 atm) . Under these conditions 
the kinetics builds up a surface population, in which O 
and CO compete for either site type at the surface. This 
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FIG. 8: (Color online) Snapshot of the steady-state surface 
population under optimum catalytic conditions at T — 600 K 
(po 2 = 1 atm, pco = 20 atm). Shown is a schematic top view 
of the simulation area of (30 x 30) surface sites, where the 
substrate bridge sites are marked by gray stripes and the cus 
sites by white stripes. Oxygen adatoms are drawn as light 
gray (red) circles and adsorbed CO molecules as dark gray 
(blue) circles. Movies of entire kMC simulations over 500 ns 
and 1 s at various pressure and temperature conditions are 
also availablei— 



strong competition is also reflected by the large fluctua- 
tions in the surface populations with time, apparent in 
Fig. ^displaying the time evolution of the system under 
exactly these gas phase conditions. At steady-state, only 
11% of all bridge sites are on average occupied by CO, 
while the remaining 89% are covered by O atoms, leav- 
ing practically no vacancies at bridge sites at these high 
pressures. Thus, only very few CO br are available to ini- 
tiate an O cus + CO br reaction. On the other hand, the 
average CO occupation at cus sites is rather high (70%), 
cf. Fig. offering much more possibilities for the 29% 
qcus p p ma tion to engage in an O cus + CO cus reaction. 
Although the low barrier O cus + CO br elementary pro- 
cess itself exhibits very suitable properties for catalysis, 
it can thus only occur too rarely in the full concert of all 
possible events to decisively affect the catalytic function- 
ality. This emphasizes the importance of the statistical 
interplay and the novel level of understanding that is pro- 
vided by the first-principles mesoscopic studies described 
here. 

The just discussed results from the kMC simulations 
show also that explicit consideration of the evolving spa- 
tial distributions and their fluctuations are a crucial as- 
pect for the high efficiency of catalysis. In a mean field 
approach the occurrence ratio of the O cus + CO cus to the 
qcus _|_ QQbr reac tion would be determined by the ratio 
of the CO cus to CO br population available to engage with 
the O cus species in a reaction, as well as the ratio of the 
two respective process rates. Even if such an approach 
would yield the same average populations at the just dis- 
cussed pressure conditions (70% CO cus vs. 11% CO br ), 
a ratio of about 1:1 would still result for the two reac- 



tions, considering the 0.1 eV difference in the two barrier 
heights. As mentioned before, the actually calculated ra- 
tio in the kMC simulations is, however, ^3:1. Visual 
inspection of the snapshot population displayed in Fig. 
[S] points immediately at the reason for this difference. 
Despite the absence of lateral interactions among the ad- 
species, they are clearly not distributed in a random ar- 
rangement, but tend to cluster into small domains, which 
extend particularly in the direction along the bridge rows 
and cus trenches. This tendency arises out of the statis- 
tical interplay of all elementary processes, but has partly 
to do with the fact that the chemical reactions can only 
attack directly neighboring O and CO pairs. Adsorbed 
O or CO inside a formed domain can therefore only be- 
come amenable to reactions, if an efficient reshuffling of 
the populations occurs, e.g. through desorption events 
creating vacancies that are then replenished with the re- 
spective other species. 

At the high pressures corresponding to the "active re- 
gion" at 600 K all lattice sites are almost always occupied, 
preventing a reshuffling due to diffusion events and sup- 
pressing particularly a reshuffling due to oxygen adsorp- 
tion, which would require two neighboring vacant sites. 
Although the impingement rate of O2 and CO from the 
gas phase is comparable at the peak activity conditions 
(Po 2 — 1 atm, pco — 20 atm), we correspondingly ob- 
serve ~ 10 3 more CO than oxygen adsorption events. 
Reshuffling at bridge sites is further inhibited due to 
the large desorption barrier of the strongly bound ma- 
jority O br species. All in all, this leads effectively to a 
much more static bridge population compared to the cus 
population, as reflected by the much smaller fluctuations 
in Fig. 0] This helps to stabilize certain local domain 
structures, which clearly differ from a random distribu- 
tion. Most notably, almost all CO br species are rapidly 
chaperoned by CO cus in the left and right neighboring 
cus sites, cf. Fig. [S] They are then no longer accessi- 
ble to the low barrier O cus -I- CO br reaction, explaining 
the much scarcer occurrence of this reaction compared to 
the ~ 1 : 1 ratio expected in a mean field picture. Since 
approaches based on macroscopic rate equations^ rely 
inherently on such a mean field approximation and are 
therefore by definition unable to properly describe the 
correlated fluctuations observed in our kMC simulations, 
it is unlikely that they could account for the correct cat- 
alytic activity in this "active region" of the surface, even 
if they were parameterized with the same first-principles 
data employed here. 

We arrive at essentially the same picture when an- 
alyzing the "active region" of the surface at T = 
350 K. It is, however, worthwhile to point out some 
specific differences. At the lower temperature, the co- 
existence region is much narrower due to the reduced 
configurational entropy, and there is correspondingly 
only a much smaller range of environmental conditions, 
which shows appreciable catalytic activity at all, cf. 
Fig. [7| Due to the significantly decelerated rates 
of the activated surface processes, the TOFs are even 
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FIG. 9: (Color online) Snapshot of the steady-state surface 
population under optimum catalytic conditions at T — 350 K 
(Po 2 = IO -10 atm, pco = 5 10 -10 atm). Shown is a schematic 
top view of the simulation area of (30 x 30) surface sites, 
where the substrate bridge sites are marked by gray stripes 
and the cus sites by white stripes. Oxygen adatoms are drawn 
as light gray (red) circles and adsorbed CO molecules as dark 
gray (blue) circles. Movies of entire kMC simulations over 
500 ns and 1 s at various pressure and temperature conditions 
are also availablei— 

at optimum performance at this temperature still or- 
ders of magnitude lower than at 600 K, actually peak- 
ing over 10 -3 site _1 s _1 (10 12 cm _2 s _1 ) only in a 
tiny (pco-,Po 2 ) region. And since this optimum region 
corresponds now to much lower pressures, the overall 
average surface population is also lower. This is il- 
lustrated in Fig. OH by a representative snapshot at 
Po 2 = 10~ 10 atm and pco = 5 10~ 10 atm, i.e. gas 
phase conditions for which we obtain a maximum TOF 
of 4.6 IO" 3 site-V 1 (4.6 10 12 cm" 2 s- 1 ). While in the 
afore discussed "active region" at 600 K there was only a 
1% vacancy concentration even at the weaker binding cus 
sites, this is now much higher, i.e. on average 30% of all 
cus sites are empty. With 50% CO cus compared to 20% 
O cus , there is still a higher CO than O cus population 
though, as was the case at the higher temperature. Simi- 
larly, O br dominates the bridge population, however, now 
with virtually 100%. Correspondingly, the lowest barrier 
O cus + CO br reaction, which would require the presence 
of CO br species, is almost completely suppressed at this 
temperature, amounting only to a mere 0.4% of the total 
TOF. 

The overall population at 350 K is therefore much 
closer to the one expected within the "constrained ther- 
modynamic equilibrium" approach, cf. Fig. [3] This 
is due to the fact that the scarce reaction events oc- 
curring at this temperature lead only to a modest re- 
duction of the O cus population. The surface structure 
and composition becomes therewith also more similar to 
some of the terminations that have hitherto been sta- 
bilized and atomically characterized in surface science 
UHV experiments. 8 i 15 i 16 i 37 i 38 i 39 i 40 i 41 i 42 The problem with 
such studies is that depending on the specific preparation 



recipe quite different surface structures can be kineti- 
cally frozen in. These need then not to have anything to 
do with the optimum disordered and dynamic adsorbate 
composition, built up by the intricate interplay of all sur- 
face processes under steady-state conditions. While cau- 
tious interpretation of corresponding experimental data 
does in principle provide a wealth of atomic-scale infor- 
mation, direct conclusions concerning the high-pressure 
steady-state reactivity are problematic. This holds in 
particular for reactivity analyses e.g. through temper- 
ature programmed reaction spectroscopy. If they are 
based on initial surface populations that are qualitatively 
different from the real steady-state ensemble, qualita- 
tively wrong conclusions can be deduce d 52 ' 53 . 



C. Uncertainties and compensation effect 

With the advent of first-principles TOF maps as those 
shown in Figs. |S]and[2| a more detailed data base is be- 
coming available for comparison with experiments, which 
will eventually advance our microscopic understanding of 
catalysis. In order to fully appreciate the quality gained 
by the underlying first-principles statistical mechanics, it 
is, however, mandatory to qualify its accuracy. Since sev- 
eral scales are bridged, uncertainties can arise at many 
levels. In the electronic regime, we have errors in the 
energetics both due to approximations in the numerical 
setup and the exchange-correlation treatment. On the 
way to the mesoscale, it is the use of transition-state the- 
ory. Finally, at the statistical mechanics level, it is the 
suitability of the employed lattice model of the surface 
and the set of elementary processes occurring on it (e.g. 
the neglect of metal precipitates, vacancies or steps). 

While it seems natural that some degree of coarse- 
graining on the way to macroscopic functionality should 
be permissible and would thereby alleviate some of the 
errors introduced in the electronic-level description, this 
must still be carefully checked for every system. In 
the present example of CO oxidation at RuO2(110) 
we are confident that the employed lattice and pro- 
cess list captures all of the essential physics of the ac- 
tive surface at steady-state discussed here. This confi- 
dence stems from extensive theoretical and experimental 
work over the last years, which established the promi- 
nent role of the two adsorption sites, and confirmed 
in particular, that although steps and domain bound- 
aries between Ru02 grains are present, their influence is 
not significant, 7 ! 11 . 17 ! 18 ! 34 ! 37 ! 38 . 39 ! 40 ! 41 ! 42 Our error analy- 
sis concentrates therefore on the individual process rates 
entering the statistical mechanical simulations. To ac- 
count for the effective uncertainty introduced by the ap- 
proximate prefactors and activation barriers, we corre- 
spondingly reran the kMC simulations for rates that were 
systematically varied by a factor of 100, which would e.g. 
correspond to a barrier variation of ±0.2 eV at 600 K. It is 
important to note, that these variations were in all cases 
done in such a way that detailed balance, cf. eq. 10, 
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was not violated. If, for example, an adsorption rate was 
lowered (reflecting e.g. an increased activation barrier), 
so was the connected desorption rate, since the particle 
would have to overcome an equally increased barrier. 

As a first observation, our results are basically not 
affected by the uncertainty in the eight diffusion pro- 
cesses. Due to the mostly quite high barriers, most of 
these events have an extremely low process rate anyway. 
Under most environmental conditions and especially in 
the catalytically "active region" their occurrence is fur- 
ther suppressed due to the almost 100% occupancy of all 
lattice sites at all times. Correspondingly, the number 
of executed diffusion events is so low, that we could in 
fact completely switch them off in the kMC simulations, 
without observing a significant variation in our results. 
Our quantitative simulations can therefore not support 
the speculations on a specific role of the O cus species pri- 
marily in replenishing once formed bridge vacancies 
Under high-pressure steady-state conditions, such diffu- 
sion events can basically not occur, highlighting again 
that an incautious extrapolation of interpretations de- 
rived from UHV experiments can be quite misleading. 

As for the remaining 18 processes, or better 9 rates plus 
detailed balance, the central result of our error analyses is 
that the central properties discussed here, e.g. the overall 
structure of the steady-state surface populations and the 
position of optimum catalytic efficiency in {T,po 2 ,pco) 
space, are surprisingly insensitive to variations in the in- 
dividual rates. This is mostly due to the fact that these 
properties are not determined by one singular process 
alone, but by the action of many players. This puts im- 
mediately more emphasis on the relative rate differences 
between the various contributing processes instead of the 
absolute values, and systematic errors in the rates like 
e.g. a general overbinding of a given exchange-correlation 
functional are then less consequential. A too high desorp- 
tion barrier goes then e.g. hand in hand with a too high 
reaction barrier, so that the reduced rate for one reaction 
process is partially compensated by an increased number 
of available surface species, leaving the net number of 
executed reactions relatively unaffected. It is important 
to realize that this implies the necessity of a consistent 
treatment of all participating processes, and a combina- 
tion of different calculations (employing different approx- 
imations) or of theoretical and experimental parameters 
could have spoiled the description. Interestingly, a recent 
study by Honkala et al. arrived at quite similar conclu- 
sions, when studying ammonia synthesis at a ruthenium 
nanoparticle with a comparable first-principles statistical 
mechanics approach^ 

We therefore find that the variations of the rates led at 
most to shifts of up to 2 orders of magnitude in transition 
pressures from one surface structure to the other, which 
does not affect the overall structure of the steady-state 
population plots in Fig. and at all. Similarly un- 
touched is the intimate connection between the narrowly 
peaked region of high catalytic activity and the realiza- 
tion of a disordered dynamic surface "phase" , where the 



kinetics builds a surface composition in which O and CO 
occupy both site types. Under these conditions the cat- 
alytic activity results from the concerted action of many 
processes, in fact always at least 9 out of our total of 26 
processes are executed with significant frequencies. This 
strong interplay makes the resulting optimum mix at the 
catalyst surface highly adaptive and renders the corre- 
sponding maximum TOFs in the core region of highest 
catalytic activity quite insensitive to modest errors in the 
rates. Even when one or more individual process rates 
were varied by a factor of 100, this led only to changes 
in the peak TOF by a factor of 4. When one reaction 
barrier was e.g. increased, either the spatial distribu- 
tion of the surface populations adapted to yield a quite 
unchanged contribution of the corresponding process to 
the total TOF, or one or several of the other reaction 
processes occurred more often and led in this way to a 
virtually unchanged net TOF. 

These results are at variance with the frequently em- 
ployed concept of one single rate-determining step, in- 
troduced initially to simplify the analysis of a reaction 
network^ Under the conditions of high catalytic ac- 
tivity, there is in the present system no high degree of 
rate control^S^ by one specific process, since the system 
has at any time several alternatives along which to pro- 
ceed. We believe that this is a quite generic feature of 
catalytic systems. Only a first-principles statistical me- 
chanics study, as the one described above, reveals these 
effects. Crucial is the consideration of a high number 
of elementary processes (compared to typical rate equa- 
tion approaches focusing on few generic representatives 
of process types) and the consideration of the explicit 
correlated spatial distribution at the surface (compared 
to the mean field approximation underlying rate equa- 
tion approaches). Although perhaps unexpected^, the 
abundance of the CO cus — O cus nearest neighbor pairs 
(as well as CO br — O cus nearest neighbor pairs) is then 
apparently playing a role of similar importance as the 
energy barriers of individual processes. In this sense, the 
frequently requested chemical accuracy for the descrip- 
tion of individual processes appears to be a misleading 
concept. At least for the present system, a careful com- 
bination of DFT and statistical mechanics is at least as 
important as good data for the energy barriers. 



IV. SUMMARY 

In summary we have presented a first-principles sta- 
tistical mechanics approach to quantitatively study 
the steady-state situation of heterogeneous catalysis. 
Density-functional theory is employed together with 
transition state theory to accurately obtain the energetics 
of all relevant elementary processes, subsequently solving 
the statistical mechanics problem by kinetic Monte Carlo 
simulations. The key merits of this two-step approach are 
that microscopic insight into the system can be gained 
by following its full dynamics from picoseconds up to 
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seconds, while simultaneously accounting for all atomic- 
scale correlations, fluctuations, and spatial distributions 
at the catalyst surface. 

In its application to the CO oxidation at the 
RuO2(110) model catalyst, turnover frequencies in un- 
precedented agreement with existing experimental data 
are achieved in a wide range of (T,po2,Pco) conditions 
ranging from UHV up to technologically relevant pres- 
sures of the order of atmospheres and up to elevated 
temperatures. The catalytic activity is narrowly peaked 
at (T,p) conditions, where the surface kinetics builds a 
disordered and dynamic adsorbate composition with O 
and CO competing for both prominent cus and bridge 
adsorption sites at the surface. The statistical analysis 
of the surface dynamics and of the various elementary 
processes in steady-state operation reveals several sur- 
prising results. Most notably, the chemical reaction with 
the most favorable energy barrier contributes only little 
to the overall CO2 production at maximum activity con- 
ditions. Although the process itself features very suitable 
properties for catalysis, it can only occur too rarely in the 
full concert of all possible events to decisively affect the 
observable macroscopic functionality. This strong inter- 



play of a larger number of elementary processes makes 
the resulting optimum mix at the catalyst surface highly 
adaptive and renders the corresponding maximum TOFs 
in the core region of highest catalytic activity quite in- 
sensitive to modest errors in the rates. This challenges 
established concepts like the idea of one rate-determining 
step or the necessity of chemical accuracy for the de- 
scription of individual processes, and clearly illustrates 
the new quality and novel insights gained by the modern 
first-principles statistical mechanics methodology devel- 
oped and applied here. 
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